pacman::p_load(sf, tmap, sfdep, tidyverse, knitr, plotly, Kendall)In-class Ex 2
Geospatial Analysis using sfdep
1 Loading of R Packages
The following packages are used for this exercise:
- sf for handling spatial data and geoprocessing
- tmap for creating thematic maps
- sfdep for creating space-time cubes and emerging hot spot analysis
- tidyverse as a universe of packages used for aspatial data transformation
- plotly for interactive charts
2 Importing data
import hunan geospatial shapefile and hunan_2012 aspatial dataframes:
hunan <- st_read(dsn = "data/geospatial",
layer = "Hunan")Reading layer `Hunan' from data source
`C:\haileycsy\ISSS624-AGA\In-class_Ex\ice2\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 88 features and 7 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 108.7831 ymin: 24.6342 xmax: 114.2544 ymax: 30.12812
Geodetic CRS: WGS 84hunan2012 <- read_csv("data/aspatial/Hunan_2012.csv")
GDPPC <- read_csv("data/aspatial/Hunan_GDPPC.csv")2.1 Joining the dataframes
Spatial features are added to the attribute dataframe as geometry column:
hunan_GDPPC<- left_join(hunan,
GDPPC,
by = "County")
glimpse(hunan_GDPPC)Rows: 1,496
Columns: 10
$ NAME_2 <chr> "Changde", "Changde", "Changde", "Changde", "Changde", "Cha…
$ ID_3 <int> 21098, 21098, 21098, 21098, 21098, 21098, 21098, 21098, 210…
$ NAME_3 <chr> "Anxiang", "Anxiang", "Anxiang", "Anxiang", "Anxiang", "Anx…
$ ENGTYPE_3 <chr> "County", "County", "County", "County", "County", "County",…
$ Shape_Leng <dbl> 1.869074, 1.869074, 1.869074, 1.869074, 1.869074, 1.869074,…
$ Shape_Area <dbl> 0.1005619, 0.1005619, 0.1005619, 0.1005619, 0.1005619, 0.10…
$ County <chr> "Anxiang", "Anxiang", "Anxiang", "Anxiang", "Anxiang", "Anx…
$ Year <dbl> 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014,…
$ GDPPC <dbl> 8184.00, 10995.00, 12670.00, 14128.00, 16763.00, 19817.00, …
$ geometry <POLYGON [°]> POLYGON ((112.0625 29.75523..., POLYGON ((112.0625 …3 Cloropleth
tmap_mode("plot")
tm_shape(hunan_GDPPC) +
tm_fill("GDPPC",
style = "quantile",
palette = "Blues",
title = "GDPPC") +
tm_borders(alpha = 0.5) +
tm_layout(main.title = "Distribution of GDP per capita by district, Hunan Province",
main.title.position = "center",
main.title.size = 1.2,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_compass(type="8star", size = 2) +
tm_scale_bar() +
tm_grid(alpha =0.2)
4 Computing Neighbors and Deriving Contiguity Weights
Neighbour Matrix and Queen’s Contiguity Spatial weights are calculated together using st_contiguity and st_weights:
wm_q <- hunan_GDPPC %>%
mutate(nb = st_contiguity(geometry),
wt = st_weights(nb,
style = "W"),
.before = 1)5 Local Measures of Autocorrelation
lisa <- wm_q %>%
mutate(
local_moran = local_moran(GDPPC,
nb,
wt,
nsim = 99),
# place new columns in front
.before = 1) %>%
# values are stored as list -- unnest as a separate column
unnest(local_moran)6 Create a time series cube
using spacetime() from sfdep package to create a spacetime cube object:
GDPPC_st <- spacetime(GDPPC,
hunan,
# define location column
.loc_col = "County",
# define time column
.time_col = "Year")6.1 Confirm if the new dataframe is a spacetime cube object
is_spacetime_cube(GDPPC_st)[1] TRUE7 Computing GI*
GDPPC_nb <- GDPPC_st %>%
# need to specify which field to calculate from using activate when using spacetime cube objects
activate(
"geometry"
) %>%
mutate(
nb = include_self(st_contiguity(geometry)),
wt = st_inverse_distance(nb,
geometry,
scale = 1,
alpha = 1
),
.before = 1
) %>%
set_nbs("nb") %>%
set_wts("wt")- activate() of dplyr package is used to activate the geometry context
- mutate() of dplyr package is used to create two new columns nb and wt.
- Then, we will activate the data context again and copy over the nb and wt columns to each time-slice using set_nbs() and set_wts()
- row order is very important so do not rearrange the observations after using set_nbs() or set_wts().
Note that this dataset now has neighbors and weights for each time-slice:
head(GDPPC_nb)# A tibble: 6 × 5
Year County GDPPC nb wt
<dbl> <chr> <dbl> <list> <list>
1 2005 Anxiang 8184 <int [6]> <dbl [6]>
2 2005 Hanshou 6560 <int [6]> <dbl [6]>
3 2005 Jinshi 9956 <int [5]> <dbl [5]>
4 2005 Li 8394 <int [5]> <dbl [5]>
5 2005 Linli 8850 <int [5]> <dbl [5]>
6 2005 Shimen 9244 <int [6]> <dbl [6]>gi_stars <- GDPPC_nb %>%
group_by(Year) %>%
mutate(gi_star = local_gstar_perm(
GDPPC, nb, wt)) %>%
tidyr::unnest(gi_star)With these Gi* measures we can then evaluate each location for a trend using the Mann-Kendall test. The code chunk below uses Changsha county.
cbg <- gi_stars %>%
ungroup() %>%
filter(County == "Changsha") |>
select(County, Year, gi_star)cbg %>%
summarise(mk = list(
unclass(
Kendall::MannKendall(gi_star)))) %>%
tidyr::unnest_wider(mk)# A tibble: 1 × 5
tau sl S D varS
<dbl> <dbl> <dbl> <dbl> <dbl>
1 0.485 0.00742 66 136. 589.In the above result, sl is the p-value. This result tells us that there is a slight upward but insignificant trend.
p <- ggplot(data = cbg,
aes(x = Year,
y = gi_star)) +
labs(title = "Changsha GI* By Year") +
geom_line() +
theme_light()
ggplotly(p)8 Emerging Hotspot analysis
ehsa <- gi_stars %>%
group_by(County) %>%
summarise(mk = list(
unclass(
Kendall::MannKendall(gi_star)))) %>%
tidyr::unnest_wider(mk)emerging <- ehsa %>%
arrange(sl, abs(tau)) %>%
slice(1:5)Lastly, we will perform EHSA analysis by using emerging_hotspot_analysis() of sfdep package. It takes a spacetime object x (i.e. GDPPC_st), and the quoted name of the variable of interest (i.e. GDPPC) for .var argument.
The k argument is used to specify the number of time lags which is set to 1 by default. Lastly, nsim map numbers of simulation to be performed.
ehsa <- emerging_hotspot_analysis(
x = GDPPC_st,
.var = "GDPPC",
k = 1,
nsim = 99
)What is the distribution of EHSA classes?
ggplot(data = ehsa,
aes(x = classification)) +
geom_bar()
9 Geographic visualisation of EHSA
hunan_ehsa <- hunan %>%
left_join(ehsa,
by = join_by(County == location))ehsa_sig <- hunan_ehsa %>%
filter(p_value < 0.05)
tmap_mode("plot")
tm_shape(hunan_ehsa) +
tm_polygons() +
tm_borders(alpha = 0.5) +
tm_shape(ehsa_sig) +
tm_fill("classification") +
tm_borders(alpha = 0.4)